---
title: "Exploratory Factor Analysis"
author: "David Barron"
date: "Trinity Term 2018"
fontsize: 10pt
output: 
  beamer_presentation:
    theme: "Madrid"
    toc: false
    slide_level: 3
    keep_tex: false
    fig_caption: false
    df_print: kable
 #   includes:
#      in_header: nup.tex
    
---


```{r setup, echo=FALSE, message=FALSE, warning=FALSE}
knitr::opts_knit$set(width = 100)
knitr::opts_chunk$set(comment =  '', echo=FALSE)
options('show.signif.stars'=FALSE)
options('digits' = 3)
#options(scipen=5)
library(ggplot2)
library(foreign)
library(effects)
library(dplyr)
library(psych)
```

# Introduction

### Introduction

Exploratory factor analysis (EFA), often just called "factor analysis," is intended to help researchers understand the structure of complex data, and possibly reduce the complexity by combining items that can be considered to be indicators of the same underlying concept.  

EFA refers to a set of procedures intended to determine the number of distinct constructs needed to account for the pattern of observed correlations among a set of measures.  These constructs are called **factors** (or more correctly "common factors"), or sometimes latent variables as distinct from the measured variables that are actually observed.

The basic premise is that observed correlations among measured variables are the result of the fact that they are influenced by the same common factor.  Hence, if we understand the relationships between common factors and measured variables, we can potentially reduce the complexity of the observed data.

The EFA model also assumes that observed measures are influenced by __unique factors__, representing all sources of variance in each measured variable that is not due to common factors.  One such source is measurement error, and indeed this is generally how unique factors are conceptualised, as if there are systemmatic, unique factors that would suggest a problem with the measure.

### Partitioning variance

$$
\text{observed variance} = \text{common variance} + \text{unique variance}
$$
$$
\text{unique variance} = \text{specific variance} + \text{error variance}
$$
An important way of assessing an EFA model is the porportion of the total variance accounted for by the common factors, which is called the **communality.**
$$
\text{communality} = \frac{\text{common variance}}{\text{observed variance}}
$$

# Basic common factor model

### Two common factors

Graphical representation of a common factor model with two orthogonal common factors and six measured variables.
![Two Common Factors](BasicEFAPath.pdf)

### Path diagramme conventions

It is conventional to represent factors/latent variables by circles or ovals, measured variables by squares or rectangles, linear causal relations by single-headed arrows, and covariances and variances by double-headed arrows.  Some relations and variances are constrained to have the value 1.0.  The common factors are not correlated with each other (i.e., they are _orthogonal_).  The unique factors are independent of each other.

### Correlation Structure Model

$$
P = \Lambda \Phi \Lambda^T + D_\psi
$$
$P$ is the observed correlation matrix (a $6 \times 6$ matrix in our example):
$$
P = \begin{pmatrix}
\rho_{11} \\
\rho_{21} & \rho_{22} \\
\rho_{31} & \rho_{32} & \rho_{33} \\
\rho_{41} & \rho_{42} & \rho_{43} & \rho_{44} \\
\rho_{51} & \rho_{52} & \rho_{53} & \rho_{54} & \rho_{55} \\
\rho_{61} & \rho_{62} & \rho_{63} & \rho_{64} & \rho_{65} & \rho_{66} \\
\end{pmatrix}
$$
$\Lambda$ is the matrix of _factor loadings_, the strength and direction of the linear relationships between common factors and measured variables:
$$
\Lambda = \begin{pmatrix}
\Lambda_{11} & \Lambda_{12} \\
\Lambda_{21} & \Lambda_{22} \\
\Lambda_{31} & \Lambda_{32} \\
\Lambda_{41} & \Lambda_{42} \\
\Lambda_{51} & \Lambda_{52} \\
\Lambda_{61} & \Lambda_{62} \\
\end{pmatrix}
$$

### Correlation structure model, continued

$\Phi$ is the correlation matrix among the common factors:
$$
\Phi = \begin{pmatrix}
1.00 & \Phi_{12}\\
\Phi_{21} & 1.00 \\
\end{pmatrix}
$$
In the orthogonal model, $\Phi_{21} = \Phi{12} = 0$, so this matrix can be omitted.  $D_\Psi$ is the correlation matrix among the unique factors.
$$
D_\Psi = \begin{pmatrix}
D_{\Psi_{11}} & 0 & 0 & 0 & 0 & 0 \\
0 & D_{\Psi_{22}} & 0 & 0 & 0 & 0\\
0 & 0 & D_{\Psi_{33}} & 0 & 0 & 0 \\
0 & 0 & 0 & D_{\Psi_{44}} & 0 & 0 \\
0 & 0 & 0 & 0 & D_{\Psi_{55}} & 0 \\
0 & 0 & 0 & 0 & 0 & D_{\Psi_{66}} \\
\end{pmatrix}
$$

### Implications

Consider two specific cases:
$$
\rho_{11} = 1.0 = \Lambda_{11}\Lambda_{11} + \Lambda_{12}\Lambda_{12} + D_{\Psi_{11}}.
$$
This implies that the variance of measured variable 1 comes from three sources, and the first term in the sum corresponds to the proportion of the variance due to common factor 1, the second term to the proportion of the variance due to common factor 2 and the third term to that due to the unique factor.

$$
\rho_{21} = \Lambda_{21}\Lambda_{11} + \Lambda_{22}\Lambda_{12}.
$$
This implies that the correlation between measured variables 1 and 2 is the sum of the product of common factor 1's loadings on MV1 and MV2, and the product of common factor 2's loadings on MV1 and MV2.  You can see that if both MV1 and MV2 load highly on either common factor 1 or 2 (or both), then they will have a high correlation.

### Data considerations

- When designing research with the expectation that EFA will be used, it is generally a good idea to plan to have at least 5 MVs per construct.  
- Keep measurement error to the minimum.
- MVs should be interval variables (or a close approximation).
- Sample sizes depend on the properties of the data. 
    - Under optimal conditions (low measurement error, 3--5 MVs per construct), sample sizes of 100 should suffice.
    - Under moderately good conditions, a sample size of 200 is adequate.
    
### Common Factor Model or Principal Components Model?

\small
There is a lot of confusion between these two.  In fact, many people report doing "factor analysis" when they in fact use principal components analysis (PCA).  It is widely believed that PCA is a "type" of factor analysis, but this is not the case.  Although it may be true that PCA often produces similar results to CFA, that is certainly not always true.  The main differences are:

- PCA was originally intended to reduce a set of MVs to the smallest set of scores that preserve as much information as possible.  That is, it is about modelling the variances in MVs, not covariance among them.
- Principal components were not intended to be thought of as corresponding to meaningful latent constructs.
- Principal components do not distinguish between common and unique sources of variance in MVs.

You can think of PCA as being represented by:
$$
P = \Lambda\Lambda^T, 
$$
that is, assuming that there is no variance due to unique factors.

The attraction of PCA has been ease of computation, but this is no longer a serious concern.  I generally prefer CFA primarily because it maps on to the idea that MVs are associated with theoretical, unobserved constructs.  

### Example data

\scriptsize
How much of each of nine types of social support do each of 200 subjects, who are graduate students, provide to friends and family in an average week?  (1) the number of hugs the person gives (Hugs), (2) the number of compliments the person gives (Comps), (3) the  number of times the person gives another person advice about  his or her personal life (PerAd), (4) the number of times the person invites someone  to social activities (SocAc), (5) the  number of times the person provides some type of professional 
advice (ProAd), (6) how often the person participates in  communal study sessions (i.e., studying in a group; ComSt),  (7) the number of times the person provides some form of  physical help, such as garden or house maintenance (PhyHlp), (8)  how often the person explicitly encourages others (Encour), and (9) how often the person tutors other students on an academic subject (Tutor).
```{r}
lower <- '1.000,
          0.666, 1.000,
          0.150, 0.247, 1.000,
          0.617, 0.576, 0.222, 1.000,
          0.541, 0.510, 0.081, 0.409, 1.000,
          0.653, 0.642, 0.164, 0.560, 0.667, 1.000,
          0.473, 0.425, 0.091, 0.338, 0.734, 0.596, 1.000,
          0.549, 0.544, 0.181, 0.448, 0.465, 0.540, 0.432, 1.000,
          0.566, 0.488, 0.120, 0.349, 0.754, 0.672, 0.718, 0.412, 1.000'
efa.eg <- lavaan::getCov(lower, names = c('Hugs', 'Comps', 'PerAd', 'SocAc', 'ProAd', 'ComSt', 'PhyHlp', 'Encour', 'Tutor'))
corrr::fashion(efa.eg, leading_zeros = TRUE)
```


# Number of common factors

### How many common factors?

Although there are likely to be theoretical expectations about the number of common factors, it is important that this be checked rather than assumed.  To obtain estimates of factor loadings in practice we have to specify in advance how many common factors to include in the model, and there are a range of methods intended to assess the appropriateness of a choice.  When making a choice, we need to find a model that:

1. Does a good job of accounting for correlations among MVs;
2. Would do substantially worse with one fewer common factor;
3. Would not do substantially better with one more common factor;
4. All common factors can be readily interpreted and related to theoretical constructs.

### Eigenvalues greater than one

A very common method is based on the eigenvalues of the sample correlation matrix.  The rule (sometimes called the _Kaiser criterion_) is that the number of common factors required is equal to the number of eigenvalues greater than 1.  This rule has some logic when it comes to PCA, but less so to CFA.  It also performs poorly in simulation studies.

The eigenvalues of the sample correlation matrix, $\theta$, are given by:

$$
P = V \theta V^{-1}
$$

### Scree test

This involves producing a graph plotting the eigenvalues in descending order.  The graph is examined to find the number of eigenvalues before the last major drop.  This gives the number of factors to be selected.  In this example, the last major drop follows the second largest eigenvalue, and so a two factor model would be used.  It seems ad hoc, but in practice often performs quite well.

```{r, out.width="65%"}
scree(efa.eg, pc = FALSE, main = '')
efa0 <- fa(efa.eg, 1, n.obs = 200, rotate = 'varimax', fm = 'ml')
efa1 <- fa(efa.eg, 2, n.obs = 200, rotate = 'varimax', fm = 'ml')
efa2 <- fa(efa.eg, 3, n.obs = 200, rotate = 'varimax', fm = 'ml')
```

### Parallel analysis

This method compares eigenvalues obtained from the sample correlation matrix with eigenvalues obtained from random data.  Random data is generated with the same sample size and number of MVs and this is used to obtain a correlation matrix.  The random data eigenvalues are typically averages over a number of simulations.  The number of eigenvalues from the sample data that are larger than the corresponding eigenvalues from random data is the number of common factors.
```{r, fig.keep = 'none', results = 'hide', message=FALSE, warning=FALSE}
fap <- fa.parallel(efa.eg, fa = 'fa', n.obs = 200, main = '')
```
```{r, out.width="65%"}
fap.dta <- data.frame(Sample = fap$fa.values, Simulated = fap$fa.sim, ID = 1:9)
fap.dta <- reshape2::melt(fap.dta, id.var = 'ID')
ggplot(fap.dta, aes(x = ID, y = value, colour = variable, shape = variable)) + theme_bw() +
  geom_line() + geom_point() + ylab('Eingenvalues') + 
  xlab('Factor number') + scale_x_continuous(breaks = 1:9) + 
  theme(legend.title = element_blank(), legend.position = c(0.6, 0.7))
```

### Likelihood ratio test

If we use maximum likelihood methods to obtain estimates of the factor loadings, then it is possible to compare the goodness of fit of models with different numbers of common factors using a likelihood ratio test.

**One-factor v. two-factor**
\tiny
```{r}
anova(efa1, efa0)
```
\normalsize

**Two-factor v. three-factor**
\tiny
```{r}
anova(efa2, efa1)
```

### RMSEA

Using a model fit index (more widely associated with the SEMs we will discuss next week).  The most common is the RMSEA (root mean squared error of approximation):

Number of factors  |  RMSEA
-------------------|---------------
1                  | `r efa0$RMSEA[1]`
2                  | `r efa1$RMSEA[1]`
3                  | `r efa2$RMSEA[1]`

A general rule of thumb is that RMSEA lower than 0.08 are acceptable, and ideally it should be below 0.05.

#### RMSEA definition

$\sqrt{[\chi^2_{model} - df_{model}]/[(N - 1) df_{model}]}$



# Rotation

### Rotation

Solutions to the factor analysis equations (i.e., sets of factor loadings) are not unique.  Although the solutions found by the fitting algorithms are the best fit, there are an infite number of solutions that have the same fit.  So, how to select the best one?  The process of finding the best set of loadings is generally called _rotation_ because you can think of it as rotating the axes that represent the common factors about the origin.  The "best" rotation will meet these criteria:

- Each row should contain at least one value very close to 0;
- Each column should contain at least $m$ values close to 0, where $m$ is the number of common factors;
- Every pair of columns should have several rows with near zero values in one column but not the other;
- Every pair of columns should have only a few rows with non-zero values in both columns.



### Plot of factor loadings

```{r, out.width="90%"}
efa1.unrot <- fa(efa.eg, 2, n.obs = 200, rotate = 'none', fm = 'ml')
efa1.ob <- fa(efa.eg, 2, n.obs = 200, rotate = 'oblimin', fm = 'ml', normalize = TRUE)
fa.dta <- data.frame(cbind(loadings(efa1.unrot), loadings(efa1), loadings(efa1.ob)))
names(fa.dta) <- c('Unrotated.1', 'Unrotated.2', 'Varimax.1', 'Varimax.2', 'Oblimin.1', 'Oblimin.2')
fa.dta$MV <- 1:9
fa.dta <- reshape2::melt(fa.dta, id.var = 'MV')
fa.dta <- tidyr::separate(fa.dta, variable, c('Rotation', 'Common.factor'))
fa.dta <- reshape2::dcast(fa.dta, ... ~ Common.factor)
names(fa.dta) <- c('MV', 'Rotation', 'CF1', 'CF2')

ggplot(fa.dta, aes(x = CF1, y = CF2, colour = Rotation, shape = Rotation)) + #geom_point() +
  theme_bw() + geom_hline(yintercept = 0) + geom_vline(xintercept = 0) + 
  geom_text(aes(label = MV)) + coord_equal()

```



### Factor loadings

\tiny
```{r}
efa1.ob
```

# Example

### Maslach burnout inventory example

```{r}
maslach <- read.dta('C:\\Users\\dbarron\\Dropbox\\Nurses survey\\maslach.dta', convert.factors = FALSE)

maslach2 <- dplyr::select(maslach, num_range('l', 1:3, width = 1), l4r, l5, l6, l7r, l8, l9r, l10, 
                          l11, l12r, num_range('l', 13:16, width = 2), l17r, l18r, l19r, l20, l21r, l22)

maslach2 <- maslach2[,c('l1','l2','l3','l8', 'l12r', 'l13', 'l14', 'l20',
                        'l5', 'l6', 'l10', 'l11', 'l15', 'l16', 'l22',
                        'l4r', 'l7r', 'l9r', 'l17r','l18r', 'l19r', 'l21r')]
maslach.r <- cov(maslach2, use = 'complete.obs')
fa.parallel(maslach.r, fa = 'fa', n.obs = 2899)
maslach.all <- fa(maslach.r, n.obs = 2899, fm = 'ml', rotate = 'oblimin', nfactors = 3)

```

### Factor loadings

\tiny

```{r}
loadings(maslach.all)
```

RMSEA: `r maslach.all$RMSEA[1]`

### Factor loadings plot

```{r}
plot(maslach.all)
```
